function verify_me
close all
% problem parameters
delta = 1e-2;
u0 = [0; 1];

T = 1;
N = 500;

% A = @(tau) [
%     2 * exp(-tau) - 0.1 -1;
%     0 -2 * exp(-tau) - 0.1];
A = @(tau) [    
    -tau.^3 ./ (1 + tau.^3), tau.^2 ./ (1 + tau.^3);
    tau ./ (1 + tau.^3), -1 ./ (1 + tau.^3)
    ];

tau = linspace(0, T, N);
dim = size(A(0),1);

% compute numerical soln
opts = odeset('abstol',1e-8, 'reltol',1e-8);
num = ode15s(@(t, X) A(delta * t) * X, [0, T / delta], u0, opts);
fprintf('num soln computed\n');
num_i = deval(num, tau / delta);


% loop over tau to get eigenvalues/eigenvectors
for i = 1:N
    
    [evec_tmp, ev_tmp] = eigs(A(tau(i)));
    
    v(:,:,i) = fliplr(evec_tmp);
    lambda(:,i) = flipud(diag(ev_tmp));
    adj(:,:,i) = inv(v(:,:,i)');
    
    c_num(:,i) = adj(:,:,i)' * num_i(:,i);
    
end

% interpolate eigenvectors for easy differentiation
v_i = @(t_i) interp3(1:dim, 1:dim, tau, v, (1:dim)', (1:dim), t_i);

% compute the c_i
c_i = adj(:,:,1)' * u0;

% loop again to compute gamma
for i = 1:N
   
    if (i == 1)
        dv = (v_i(tau(i) + 1e-5) - v_i(tau(i))) / 1e-5;
    elseif (i == N)
        dv = (v_i(tau(i)) - v_i(tau(i) - 1e-6)) / 1e-6;
    else
        dv = (v_i(tau(i) + 1e-5) - v_i(tau(i) - 1e-5)) / 2e-5;
    end
    
    % gamma(j,k,.) = v*(j)' * v(k)
    gamma(:,:,i) = -adj(:,:,i)'*dv;
    
end
% g_i = @(t_i) interp3(1:dim, 1:dim, tau, gamma, (1:dim)', (1:dim), t_i);
% g_i(1/2)

% compute beta
tmp = zeros(N, dim);
for i = 1:dim
    tmp(:,i) = reshape(gamma(i,i,:), N, 1);
end
beta = cumtrapz(tau', lambda' / delta + tmp)';

% compute leading order coeffs
c_lo(1,:) = c_i(1) * exp(beta(1,:));
c_lo(2,:) = c_i(2) * exp(beta(2,:));


% do the corrections

int = cumtrapz(tau, reshape(gamma(1,2,:) .* gamma(2, 1,:), 1, N) ./ (lambda(1,:) - lambda(2,:)));

c_corr(1,:) =-c_i(2) * reshape(gamma(1,2,:), 1, N) ./ (lambda(1,:) - lambda(2,:)) .* exp(beta(2,:)) + ...
    (c_i(1) * int + gamma(1,2,1) * c_i(2) / (lambda(1,1) - lambda(2,1))) .* exp(beta(1,:));

c_corr(2,:) = c_i(1) * reshape(gamma(2,1,:), 1, N) ./ (lambda(1,:) - lambda(2,:)) .* exp(beta(1,:)) - ...
    (c_i(2) * int + gamma(2,1,1) * c_i(1) / (lambda(1,1) - lambda(2,1))) .* exp(beta(2,:));



% plot the coefficients
subplot(2,1,1);
plot(tau, c_num(1,:), 'k', tau, c_lo(1,:), 'b', tau, c_lo(1,:) + delta * c_corr(1,:), 'r');
subplot(2,1,2);
plot(tau, c_num(2,:), 'k', tau, c_lo(2,:), 'b', tau, c_lo(2,:) + delta * c_corr(2,:), 'r');


% plot the relative error
figure
subplot(2,1,1);
semilogy(tau, abs((c_lo(1,:) - c_num(1,:)) ./ c_num(1,:)) , 'b', tau, abs((c_lo(1,:) + delta * c_corr(1,:) - c_num(1,:)) ./ c_num(1,:)), 'r');
hold on;
semilogy(tau, ones(1,N) * delta, 'k--');
semilogy(tau, ones(1,N) * delta^2, 'k--');
subplot(2,1,2);
semilogy(tau, abs((c_lo(2,:) - c_num(2,:)) ./ c_num(2,:)) , 'b', tau, abs((c_lo(2,:) + delta * c_corr(2,:) - c_num(2,:)) ./ c_num(2,:)), 'r');
hold on;
semilogy(tau, ones(1,N) * delta, 'k--');
semilogy(tau, ones(1,N) * delta^2, 'k--');